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We show by unbiased numerical calculations that the ferromagnetic nearest-neighbor exchange 
interaction stabilizes a vector spin chiral order against the quantum fluctuation in a frustrated spin-| 
chain relevant to multiferroic cuprates, LiCu2 02 and LiCuV04. Our realistic semi-classical analyses 
for LiCu2 02 resolve controversies on the helical magnetic structure and unveil the pseudo-Nambu- 
Goldstone modes as the origin of experimentally observed electromagnons. 
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The parity symmetry may be broken spontaneously in 
magnets [1-3] even if the original crystal structure pre- 
serves it. It occurs, for instance, when the vector chi- 
rality [3] Sr') of neighbouring spins Sr 

and Sr' acquires a nonzero macroscopic average in ge- 
ometrically frustrated magnets [Fig. 1(a)]. The helical 
magnetism [1, 2] in which the spins align in a helix with 
a particular handedness is a typical example. This issue 
of the spin chiral order [3] and the associated electromag- 
netic excitations have been highlighted by a recently dis- 
covered multiferroic behavior, a ferroelectricity induced 
by a spin cycloid [4, 5]. In spite of theoretical advances 
in the understanding of the static magnetoelectric effect 
[6-8], the dynamical effects [9-12] remain controversial. 

In helical magnets, three Nambu-Goldstone modes ap- 
pear in principle [1, 2] since the SU(2) symmetry of the 
spins are fully broken. In the case of the spin cycloid 
[Fig. 1(b)], one corresponds to a phason describing an in- 
finitesimal change of the pitch of the cycloid [Fig. 1(c)]. 
The other two represent infinitesimal rotations of the cy- 
cloid plane [Fig. 1(d)]. Because of magnetic anisotropy, 
they usually acquire an energy gap and can then be 
probed with the antiferromagnetic resonance [13]. These 
modes may also be excited by the electric component of 
light through the magnetoelectric coupling. Recent THz 
spectroscopy experiments on i?Mn03 [9, 10] have demon- 
strated that this electromagnon spectrum grows below 
the ferroelectric transition temperature. However, be- 
cause of the large GdFeOa-type distortion, the observed 
spectrum is dominated by high-energy magnons at the 
Brillouin-zone boundary [10] through the magnetostric- 
tion mechanism [14]. The roles of the Nambu-Goldstone 
modes remain open [11]. 

The ferroelectricity associated with the spin cycloid 
has also been found in quasi-one-dimensional (QID) 
cuprates, LiCu202 [15] and LiCuV04 [16]. These spin- 
chain compounds are advantageous in observing the 
lowest-energy magnons in the optical spectrum since the 
crystal structure is less distorted along the chains. In 
particular, in LiCu202, the contribution from the zone- 
boundary magnons along the chain is prohibited by the 
symmetry. However, the understanding of experimen- 
tal findings on the QID multiferroic cuprates remains 



controversial, including the magnetic ordering pattern 
[12, 17-20] and the dynamical properties [12] in LiCu202. 
Furthermore, it is by far nontrivial from the theoretical 
viewpoint whether the chiral LRO parasitic to the heli- 
cal magnetism can be stabilized against quantum fluctu- 
ations. In fact, strong quantum fluctuations in one di- 
mension often lead to a valence-bond order accompanied 
by a moderately large gap in the spin excitations [21, 22]. 
This spin gap can prevail over weak inter-chain couplings 
and prevent the helical magnetic long-range order (LRO), 
as in CuGeOa [23]. A possibility that the chiral ordered 
phase can appear for a weak easy-plane anisotropy in the 
frustrated spin-^ chain has been addressed [24, 25], and 
conclusive calculations are called for. 

In this Letter, we develop a comprehensive theory for 




FIG. 1: (Color) (a) A ID frustrated spin model with Ji and 
J2{> 0) being the nearest-neighbor and second-neighbor ex- 
change couplings, (b) A spin cycloid in the propagation direc- 
tion of the J1-J2 spin chain. Spins rotate within the common 
plane represented by the disks. (c,d) Three Nambu-Goldstone 
modes in the SU(2)-symmetric model; a phason [(c)] and two 
rotation modes of the spiral plane [(d)]. (e) Cu^'*', Cu2^'*' and 
ions in LiCu2 02, and inter-chain Heisenberg exchange 
{Ja,±) and inter-layer Dzyaloshinskii-Moriya interactions. On 
the arrows with four different colors, four distinct DM vec- 
tors {D") are assigned, (f) Hypothetical magnetic structure 
(red arrows) and lowest-energy spin-excitation mode (rota- 
tions around blue arrows) in the THz range for LiCu2 02. 
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FIG. 2: (Color) Ground-state phase diagram of the spin chain model Hid- (a) The classical phase diagram and relevant QID 
spin- 1 materials. Estimates of J1/J2 for materials are taken from Refs. 17, 26, 27. The materials shown in green exhibit an 
antiferromagnetic LRO, but their detailed magnetic structures have not been established yet. In Rb2Cu2Mo30i2, a magnetic 
LRO has not been detected down to 2 K [26]. Usually, 1 — A < 0.05 in cuprates. (b) The quantum phase diagram and the map 
of the chiral order parameter for S — ^. The horizontal axis J1/J2 is the same as in the panel (a). The boundary between 
the dimer phase and the chiral phase in the case of Ji > agrees with the previous numerical study [22]. For IJ1I/J2 > 4, 
Tomonaga-Luttinger liquid (TLL) phases appear. Detailed methods of identifying each phase and the phase boundaries are 
explained in Supplementary Material. It was difficult to determine the phase boundary for small | J1I/J2 (dashed line). 



QID multiferroic cuprates. We uncover that the ferro- 
magnetic nearest-neighbor exchange coupling Ji stabi- 
lizes a chiral order in the frustrated spin-^ chain un- 
der weak easy-plane anisotropy, in sharp contrast to 
the case of the antiferromagnetic Ji [Fig. 2(b)]. With 
weak three-dimensional (3D) couplings, this roughly con- 
trols whether the ground state exhibits the helical mag- 
netic order or the Neel/dimer order. Our phase diagram 
is useful for classifying several QID cuprates [26, 27] 
[Fig. 2] . We also theoretically clarify the magnetic order- 
ing structure and the electromagnetic excitation spectra 
in LiCu2 02, which show overall agreements with exper- 
iments [12, 15, 17-20]. These agreements give evidence 
that weak inter-chain couplings and magnetic anisotropy 
allow the electric component of light to excite otherwise 
zero-energy Nambu-Goldstone modes. 

We first reveal the origin of the vector chiral order un- 
derlying the spin spiral in LiCuV04 and LiCu202. At 
higher temperatures than the inter-chain couplings of or- 
der of the Neel temperature (2.5 K and 24 K for LiCuV04 
and LiCu202, respectively), they are described in terms 
of a one-dimensional (ID) spin-^ model [Fig. 1(a)], 

Hid = ^ Jn ^ {S^S''!I_^^^ -f SjS^^,^ + l^S^S^^^ , 

n=l,2 i 

(1) 

( J2 > 0) with the electronic spin S-j = (Sj, S'J, S^) at the 
Cu^"*" site j in the chain, and the small symmetric easy- 
plane exchange anisotropy 0<1 — A^l. In the classi- 



cal limit, a helical magnetic ground state is realized for 
\Ji\/J2 < 4 [Fig. 2(a)], with a finite uniform vector spin 
chirality {kj = . In the quantum spin-^ case, this 
chiral order tends to yield to valence bond orders [21, 22]. 
However, our previous finite-size calculation of the spin 
Drude weight, supplemented by the bosonization analy- 
sis, has suggested that the chiral phase might survive in 
a wider range for ferromagnetic Ji [25]. To precisely de- 
termine the ground-state phase diagram in a conclusive 
manner, we have performed numerical calculations based 
on the time evolving block decimation algorithm for an 
infinite system (iTEBD) [28]. Figure 2(b) presents the 
global phase diagram and a profile of the chiral order pa- 
rameter in the space of J1/J2 and A. It confirms that 
the chiral phase extends over a wide region for ferromag- 
netic Ji < 0. In particular, the chiral order is pronounced 
by small magnetic anisotropy 1 — A > and a moderate 
value of I J1I/J2 that roughly correspond to the multifer- 
roic cuprates. This chiral LRO is stable up to a close 
vicinity of the S'C/(2)-symmetric case A = 1, where it is 
replaced by a dimer order. It is also replaced with an- 
other dimer order with a unit of Wl) + by stronger 
anisotropy < A < 0.6. On the other hand, for antifer- 
romagnetic Ji > 0, the singlet-dimer order accompanied 
by the spin gap is robust [21, 22]. 

The above results clearly explain why the QID helical 
magnets or the candidate materials are found dominantly 
on the ferromagnetic side of Ji while the materials lo- 
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cated on the antifcrromagnctic side of Ji usually show a 
spin-singlet dimer order or a coUinear antiferromagnetic 
LRO. In real materials having weak inter-chain couplings, 
the helical spin quasi-LRO in the chiral phase [24] readily 
evolves into a genuine helical LRO. Even when the mag- 
netic anisotropy is so small that a material with Ji < is 
located in a narrow dimer phase near A = 1, the spin gap 
is orders of magnitude smaller than for Ji > 0. There- 
fore, the weak 3D couplings may close the spin gap and 
recover the helical magnetic LRO. Once the helical mag- 
netic LRO readily occurs at a finite temperature on the 
ferromagnetic side of Ji, as in the case of LiCu202, then 
the static and dynamical properties can be calculated 
through semi-classical analyses [2]. In the rest of this 
Letter, constructing a realistic spin model, we will give 
a theoretical explanation of overall experimental findings 
on LiCu202. 

The material contains four J1-J2 spin chains [Fig. 1(a)], 

each of which lies in one of four stacked ab layers (I to 
IV) in the unit cell [Fig. 1(e)]. Neutron-scattering exper- 
iments [17] have unveiled a magnetic Bragg peak at an 
incommensurate wavcvcctor (0.5,0.173,0) in the recipro- 
cal lattice unit, and the QID spin-wave dispersion along 
the b axis. The observation of the ferroelectric polariza- 
tion P \\ c [15, 19] and the polarized neutron-scattering 
experiment [19] indicate the emergence of the uniform be 
spin-spiral component in the lowest-temperature phase. 
Then, a minimal spin model should include the interac- 
tions shown in Fig. 1(e). In each layer, the adjacent spin 
chains along the a axis are coupled via the Heisenberg 
exchange interaction with moderate coupling strength Ja 
[Fig. 1(e)]. Two pairs of the adjacent layers, II and III, 
and IV and I, are also coupled via the nearest-neighbor 
inter-chain Heisenberg exchange coupling Jl in a zigzag 
manner [Fig. 1(e)]. On the other hand, a small Heisen- 
berg exchange coupling J' (not shown) between the other 
two pairs of layers, I and II, and HI and IV, induces 
an incommensurability along the a axis, i.e., Qa ^ 7r/a, 
which disagrees with the experiment [17, 19, 20]. To sta- 
bilize the experimentally reported uniform be spiral com- 
ponent [12, 15, 19, 20] and to pin the commensurability 
along the a axis, the a component of the Dzyaloshinskii- 
Moriya (DM) vectors for the inter-layer DM interactions 
between I and II and between HI and IV is indispens- 
able. This can be provided through the second-neighbor 
bonds characterized by the four DM vectors D" , which 
are specified by the amplitudes of the a, b, and c com- 
ponents D", D", and D" [Fig. 1(e)] from the symmetry 
argument. For simplicity, wc take J' = D" = D" = 0. 
The intra-chain DM vectors specified by the amplitude 
D and the angle 6 are uniform within each layer, but 
the symmetries of the crystal require that their direc- 
tions alternate in the layer index; Di = D{cos6,O,sm0), 
Dii = D{-cos9,0,sme), Dm = D{cose,0,-smO), 
and DiY = D{— cos9,0, — sm0). The classical anal- 
ysis of this model Hamiltonian shows that the intra- 



chain DM interaction rotates the cycloid planes about 
the b axis from the be plane by different angles ipi, (p2, 
—ifi, and —(p2 for the layers I, II, III, and IV, respec- 
tively. With the Bragg-peak position (0.5,0.173,0) [17], 
a choice of (<pi, iy92)/(27r) = (—0.08, —0.02) reproduce the 
'^Li nuclear magnetic resonance spectra [18, 20] (Sup- 
plementary Figure S4). Then, the exchange and DM 
coupling parameters are determined to reproduce the 
ordering wavevector Q, the rotation angles (i^i,</'2) of 
the spiral planes, the anisotropy in the antiferromag- 
netic resonance spectra [12]; {Ji, J2, Ja, J±, D, D'^) = 
(-11.3,5.9,1.1,0.08,0.37,0.26) meV and 6'/(27r) = 0.43. 
It is noteworthy that our choice (Ji, J2, Ja, J±) is within 
a reasonable range with respect to previous analyses of 
inelastic neutron-scattering experiments [17]. The ob- 
tained magnetic structure is presented in Fig. 1(f). 

The associated QID spin- wave dispersions are shown 
in Fig. 3(a). Magnetic anisotropy has produced the 
energy gap in otherwise zero-energy Nambu-Goldstone 
modes at q = ±Q. Since four spins exist in the unit 
cell, there appear an acoustic branch and three opti- 
cal branches. The splitting of the four branches is pro- 
nounced at the dispersion minima, q = and q = ±Q, 
where we further take account of the bilinear mixing 
of q = and q = ztQ magnons through the degener- 
ate perturbation theory of inter-layer and DM couplings. 
Figure 3(b) shows the anisotropy of the obtained anti- 
ferromagnetic resonance spectra Im[x'^°(a;) — (x""(cj) -|- 
X^''{ijj))/'2], which remarkably agrees with the experimen- 
tal result [12]. The lowest-energy mode at w ~ 0.8 meV 
is associated with a linear combination of the intra-layer 
phason modes where the change of the pitch of the cy- 
cloid in the layers I and II is opposite to that in the layers 
HI and IV [Fig. 1(f)]. It is also accompanied by a small 
fraction of 2Q oscillation. The acoustic mode was ob- 
served at ~ 30 GHz ~ 1 cm-i ~ 0.12 meV with the 
electron spin resonance experiment [13], which is much 
smaller than the THz frequency of our interest. 

Our analyses of magnetoelectric couplings (based on 
Rcf. 14) show that the dominant contribution to the di- 
electric absorption spectra arises from the vector chiral- 
ity on the bonds connecting the adjacent layers (Sup- 
plementary Material). The dominance of the chiral- 
ity contribution over the magnetostriction sharply con- 
trasts with the case of iJMnOa [10], and is ascribed to 
the fact that the unit cell contains only a single spin 
in the direction of the ordering wavevector. Through 
this magnetoelectric coupling, most of the spin-wave 
modes are clcctric-dipole active. Since main peaks ob- 
served in the THz spectroscopy [12] are much sharper 
than the magnon bandwidth of order of jJij and J2, 
wc take account only of the onc-magnon contributions. 
The anisotropy of thus determined dielectric functions, 
Im[£^^(a;) - (e««(a;) -|- s'>''{u))/2], shows a reasonable 
agreement with the experiment [12] [Fig. 3(c)]. As ex- 
pected, the lowest-energy mode described in Fig. 1(f) has 
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FIG. 3: (Color online) Semi-classical analyses for LiCu202. 
(a) Spin-wave dispersions without taking account of the bi- 
linear coupling among q — 0, ±Q modes. This coupling only 
modifies the modes around q — and ±Q, marked with green 
and blue (dashed) circles. The dispersion along the c axis is 
quite small, since J±, D'^ < J a <Si \Ji\, J2- (b,c) Anisotropy in 
the electromagnetic absorption spectra due to the magnetic 
and electric components of light, in units of (p/^s)'^ meV~^ 
and cm~^, respectively. The spectrum consists of eleven 
gapped modes after the reconstruction of the q — 0, ±Q 
modes due to their bilinear coupling; oj = 0.812, 1.431, 1.525, 
1.717 (doubly degenerate), 2.469, 2.634, 2.977, 2.983, 2.991 
(doubly degenerate) meV. Inclusion of even higher harmonics 
such as q = ±2Q modes increases the number of the modes 
and eventually broadens the peaks. The series of (5-function 
peaks are replaced with the sum of Gaussians broadened with 
the width 0.05 meV. Points with dashed lines denote the in- 
tegrated intensity experimentally observed in the THz spec- 
troscopy [12] in unit of cm~^. We note that the tilting of the 
spiral planes crucially changes the signals of electromagnons. 



significant amplitude as well. We stress that excitations 
induced by electric and magnetic components of light ap- 
pear with comparable amplitude in this system. This 
is unique to spin-^ systems, where the magnetostriction 
contribution to the optics is absent and the magnetoelec- 
tric coupling is controlled by quite a small ratio of the 
relativistic spin-orbit coupling to the on-site Coulomb re- 
pulsion. 

Dynamical magnetoelectric effects due to the Nambu- 



Goldstone modes have been identified in the dielec- 
tric spectrum of the QID spin-^ multiferroic cuprate 
LiCu202. Further experiments on QID cuprates with 
weaker 3D coupfings as in LiCuV04 [16, 29] and 
Rb2Cu2Mo30i2 [26] might also uncover gapless dielec- 
tric spectra due to phasons and chirality solitons [25]. 
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ITEBD CALCULATION 

In this section, we will give detailed explanations on 
the way how we apply the infinite time-evolving block 
decimation (iTEBD) algorithm, which has recently been 
developed by Vidal [SI], for our precise determination of 
the phase diagram of the ID spin-^ J1-J2 XXZ model. 

The method is based on a matrix product represen- 
tation of a state, as is the density matrix renormaliza- 
tion group (DMRG) algorithm [S2] . A key feature of the 
iTEBD algorithm is that it directly treats an infinite sys- 
tem by exploiting the translational invariance, and is free 
from finite-size or boundary effects. 

A matrix product state (MPS) is an efficient and well- 
controlled variational ansatz for a ID quantum groimd 
state. We consider an infinite spin-^ chain and define 
a local basis \sj) {sj =t,4,) for each site j G Z. To ob- 
tain a matrix product representation of a given state j^*), 
we first consider its Schmidt decomposition into the left 
and right half-infinite chains around a certain bond, for 
instance, (0, 1): 

a=l 

Here, [<0] and [1[>] denote two half-infinite chains rang- 
ing over the sites with j < and 1 < j respectively, 
and Aq' are the Schmidt coefficients controlling the en- 
tanglement between the two half chains. {I'S!:^"')} and 
{l^a'^')} are the orthonormal sets of states on the left 
and right half chains, respectively. In general, an infinite 
sequence of Schmidt coefficients is required to exactly 
express a given state in an infinite system. Here, we 
introduce its approximate description by truncating the 
sequence: a finite number x of the largest Schmidt coef- 
ficients are taken into account. A similar decomposition 
around the adjacent bond (1,2) gives: 

/3=1 

Now one can consider the expansion of the state in [1>] 
using a local state at the site 1 and the state in [2[>]: 

l^^=) = Er™.rf'l«i''')l<')- (S3) 



Here we have introduced a three-index tensor Pt^l. Plug- 
ging Eq. (S3) into Eq. (SI), we obtain 

I*) = E AWrW ,APl|$[<°])|.iW)|$f 1). (S4) 

a, Si ,0 

Repeating this procedure over all the bonds, the 
state is eventually represented as an infinite prod- 
uct of vectors {\a^}a=i,— ,x ^^'^ three-index tensors 

{r^.l ,;g}a,/3=i, • ,x;s,=T.i defined on the bond (j - 
and tlie site j, respectively. The Schmidt rank x controls 
the precision of the approximation, and plays a similar 
role to the number of states to be kept in the DMRG 
algorithm [S2]. With larger x, one can describe longer- 
distance correlations more correctly. A simplification is 
achieved in the presence of the translational symmetry; 
tensors X^^^ and fI-'J have a certain periodicity in j. For 
the analysis of the model Hm that includes nearest- 
neighbor and second-neighbor interactions, we employ a 
period-4 structure as in Fig. SI (a). 

The calculation goes as follows. We first prepare an ar- 
bitrary initial state j^'o) in the form of the MPS with the 
4-lattice periodicity as shown in the top line of Fig. SI (a). 
Then, we perform the imaginary-time evolution so that 
it converges to the ground state I'J'gs): 

= lim II _hJ ° II (S5) 

In the present case, the imaginary-time evolution opera- 
tor e~^^'^ for a small time interval At can be expanded 
through a Suzuki- Trotter decomposition as a sequence 
of three-site gates {U} as shown in Fig. Sl(a). We now 
consider an operation of a single gate U on three sites 
labeled by (1, 2, 3). It is useful to rewrite the state as 

C«l,Sl,S2,S3,a4 

(S6) 

with a short-hand notation on the coefficients [see 
Fig. Sl(b)], 

e = V aWfW xPirPi xi^irPi 

a2,a3 

(S7) 

After the operation of \J on the sites (1,2,3), the coeffi- 
cients are replaced by 

^aiSiS2S3a4 — ^ ] (-^l , 52 , S3 I t/|siS2S3)0a^5^s2S3a4 • 

SliS2,S3 

(S8) 
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FIG. SI: (color) Processes in the present iTEBD algorithm, (a) The period-4 matrix product state and its imaginary-time 
evolution, (b) Tensor Q encapsulating the information of three sites (1, 2, 3). (c) Two-step singular-value decompositions (SVD) 
of O carried out after operating a single imaginary-time evolution U. 



The original form of the MPS can be obtained by 
performing singular value decompositions twice for Q 
[Fig. Sl(c)]. Here we truncate the tensors taking 
only the x largest singular values (Schmidt coefficients) 
{^a}a=i,- - ,x again. We repeat this imaginary-time evo- 
lution by taking the sites (1,2,3), (2,3,4), and (4,5,6), one 
after another [Fig. Sl(a)]. Since the same operator U is 
applied for every four sites in parallel in an infinite chain 
at each imaginary time slice, the period-4 structure of 
the state is preserved. 



DETERMINATION OF THE PHASE DIAGRAM 
OF THE ID MODEL 
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Here, we explain how we determine the phase diagram 
(Fig. 2(b) of the main body) for the ID spin-i J1-J2 XXZ 
model, Eq. (1) of the main body. 

A convergence to a vector-chiral long-range ordered 
state is allowed for by including complex elements in the 
initial MPS as a "seed" for the symmetry breaking [S3]. 
The "spontaneous" chiral ordering can then be observed 
directly by calculating the local chirality 

K^- ((Si X ^2)^). (S9) 

It turned out that is uniform and independent of the 
sites inside the chiral phase shown in Fig. 2(b) of the 
main body. The chiral order parameter calculated in 
this way for the model Huj with J1/J2 = —2 is shown in 
Fig. S2. The rapid growth of the order parameter allows 
us to determine the transition points with rather high 
accuracy. The x-dependence of the data is visible only 
in the close vicinity of the transition points. Therefore 
the calculation with x — 200 is satisfactory for most of 
the parameter space, and this rank is employed in the 
calculation shown in Fig. 2(b) of the main text. On the 
antiferromagnetic side of Ji, the phase boundary deter- 
mined in this way agrees well with the previous DMRG 
result [S5]. In the XY case A = 0, the phase diagram 



FIG. S2: The chiral order parameter for J1/J2 = —2 of 
the ID spin- 1 model, calculated by iTEBD with the ranks 
X ~ 200, 300. The vertical lines indicate estimated transition 
points. 

is symmetric for Ji < and Ji > because the sign 
of Ji can be reversed under the unitary transformation. 

The dimer phases can be partly discriminated by cal- 
culating two types of dimer order parameters: 

D723 = (S^.S^ + SfSD-iS^SI + SlSi), (SIO) 
-D123 — S1S2 — S2S^. (sii) 

An alternation of the sign of D^lg and/or DI23 along the 
spin chain indicates some sort of dimer ordering. In the 
singlet dimer phase for Ji > 0, both D^23 and 1)^23 ^~ 
nite and have the same sign. In particular, for J\/ J2 = 2 
and A > —1/2, the ground state is exactly given by the 
product of nearest-neighbor singlets (| ti) ~ I it))/^/2 
[S6], which satisfies the relation — 2I?f23 — ~l/2- 

Applying the transformation W in the XY case, one finds 
that the exact ground state for J1/J2 = —2 and A = 
is given by the dimer state whose unit is replaced by the 
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triplet type, i.e., (| ti) + I it))/v^- This state satis- 
fies the relation -I^ilg = 2DI23 = -1/2. Our iTEBD 
calculation shows that -D123 ^^'^ -^123 always show op- 
posite signs inside the dimer phase appearing in A < 0.7 
on the ferromagnetic side of Ji, indicating that such a 
"triplet" character of the wavefunction survives in this 
phase. On the other hand, in the other dimer phase 
close to the isotropic case A = 1 for Ji < 0, both D^23 
and Df23 show finite values with the same sign, as in 
the singlet dimer phase appearing for Ji > 0, though 
their amplitudes are significantly reduced, as suggested 
by field-theoretical analyses [S4, S7]. 

The Tomonaga-Luttinger liquid (TLL) phases appear- 
ing for IJ1I/J2 ^ 4 have the instability toward gapped 
antiferromagnetic (Neel) and dimer phases [S8, S9]. The 
phase transitions to these two ordered phases can be an- 
alyzed efficiently by the level spectroscopy method [SIO], 
which combines the effective sine-Gordon theory with 
exact-diagonalization calculations. The phase boundary 
between the TLL phase and the dimer phase has been 
thus determined previously both for Ji > [SIO] and for 
Ji < [Sll, SI 2]. Besides, a Neel ordered phase with 
an antiferromagnetic long-range-order in the z direction 
has been found between the TLL and chiral phases in 



0.7< A<0.9 [S12]. 



CRYSTAL STRUCTURE OF LICU2O2 AND THE 
SPIN MODEL HAMILTONIAN 

Here, wc explicitly give a mathematical expression of 
the model Hamiltonian for LiCu2 02. 

LiCu202 has the space group Pnma with the lattice 
constants a = 5.726 A,b = 2.8587 A, and c = 12.4137 A, 
an inversion center, and a mirror plane perpendicular to 
the b axis. For spin-^, magnetic anisotropy is usually 
dominated by the Dzyaloshinskii-Moriya (DM) interac- 
tion. The symmetric anisotropic exchange coupling is of 
higher-order in the relativistic spin-orbit coupling. Tak- 
ing account of the symmetry properties of the crystal, we 
can construct the Hamiltonian including the DM inter- 
actions; 

IV 

'H=Y,^^r^ +'^±+'^ac, (S12) 
m=I 

where 



nj2D 

"•TO 

Hi. 

Hac 



+D'i, [Si,i,j,e X ('S'li.i-i,^ — -Sii.j/) -l- {Siii^i-ij^e — Sui^ij^^) x -S'lv.i,^,^]'' 
r 1 " 



<5a,<56=±l 



5a,St,=±l 



(S15) 



Here, Sm,i,j,e represents the electronic spin located at the 
position rm,i,j,e of the Cu^+ ion, with 

ri,i,j,t = {{i + 0.12)a, (j + 0.25)5, {£ + 0.905)c), (S16) 

ru,i,j,e = {t + 0.62)a, (j + 0.25)6, (f. + 0.595)c), (S17) 

ruu,3,£ = (i + 0.38)a, (j + 0.75)6, {£ + 0.405)c), (S18) 

riv,»,j/ = {i- 0.12)a, {j + 0.75)6, {£ + 0.095)c). (S19) 



The Heisenberg exchange couplings Ji, J2, Ja, J±, and 
J', the intra-chain DM vectors Di^...jv specified by D 
and 9, and the inter-chain DM couplings D", D", and D" 
have been explained in the main text. A finite but rather 
weak nearest-neighbor inter-chain DM vector {0,D[,0) 
may also exist because of a slight shift of the atoms from 
the inversion-symmetric positions about the bond center. 
Most likely, however, the second-neighbor inter-chain 
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DM couplings D" ^ [Fig. 1(e) of the main body] prevail 
over D'jj, since the inversion symmetry is broken more 
significantly in this geometry and the relevant exchange 
processes include common paths. In particular, £)" is the 
only DM coupling that favors the uniform be spin-spiral 
component observed in experiments [S13, S14, S15, S16], 
and should play a central role in fixing the spiral planes. 
To reduce the number of parameters, we neglect J', Z3[,, 
D'l^, and Z?" for LiCu202. Actually, the other compo- 
nents favors alternating spiral structures as proposed for 
NaCu202 [S17]. Note that the existence of the inver- 
sion center in the bonds connecting zigzag-ladder pairs, 
namely, between the layers II and III and between the 
layers I and IV, prohibits the DM interaction on these 
bonds. We ignore the DM interaction on the bonds di- 
rected to the a axis, since it is small compared with Ja 
and does not stabilize the experimentally observed stag- 
gered spin modulation along the a axis. 



CLASSICAL ANALYSES OF MAGNETIC ORDER 

Here, we explain details of the classical analysis of the 
magnetic structure for LiCu2 02, whose result is shown 
in Fig. 1(f) of the main body. We ignore an ellipticity 
caused by a density wave of magnons, which is beyond 
the scope of a linear spin-wave theory. 

We assume the magnetic structure of {Sm,i,j.e) = 



Sn 



where n 



) is a 



unit vector whose components are given by 



"m.i = cos Iprn Sm ifm COS Om,i,j,e " siu tpm Sin Om,i,j,e, 

(S20a) 

= sin V^m '^in 'P^n COS 0„i,i,j,£ + COS l/j^ sin Om,i,j/, 

(S20b) 
(S20c) 



'm.i.j.l 



COS (pm COS 6m,i,jl- 



Here, ipm and ipm represent the successive rotation an- 
gles of the cycloid plane in the mth layer about the h axis 
from the be plane and then about the c axis, respectively. 
We have introduced the phase Om,i,ji of the spin spiral 
at the site rm.ijj in the mth chain as 

ei,,,,,i = Q-{ia,jb,£c) + (f>i, (S21) 

Ou,^,J,i = Q-{{i + l/2)a,jb, £c) + 4>ii, (S22) 

Om,^,J,l = Q-{{i + l/2)a, (j + 1/2)6, &) + ^m, (S23) 

elv,^,JJ = Q ■ {ia, {j + 1/2)6, ic) + 4>iv (S24) 

We substitute these expressions into the Hamiltonian 
Eq. (S12) with D'^ = D'^ = D'^ = 0, and minimize it 
with respect to Q, ipm, ipm, and (pm {"ni — I, - ■ ■ ,IV). 
Note that large values of Z3" and Ja > compared 
with I J' I are required for reproducing the experimen- 
tally observed ordering wavevector Q ~ {Qa,Qb,0) with 
Qa = Tr/a [S14, S15, S18, S19]. In fact, once Qa = 7r/a is 



obtained, small J' does not alter the magnetic structure 
at all. Hence, we take J' = for simplicity. 

Then, the minimization procedure gives generic condi- 
tions for the rotation angles, ipm — 0, y^i — ~fiii — fi, 
and ifji = —fiv = 'i^2i a-nd those for the initial phases 
(t^ii ~ 4>i ^ -7r/2, 0/// - 4>i ^ 7r/2, and 0/v - 0/ = tt, 
as described in the text and Fig. 1(f) of the main body. 



DETERMINATION OF THE EXCHANGE AND 
DM COUPLING PARAMETERS 

Here, we explain in detail how we have determined 
the set of parameter values {Ji, J2, Ja, J±, D, D'a) = 
(-11.3,5.9,1.1,0.08,0.37,0.26) meV, and = 0.43 x 27r 
from a fitting with experimental results. 

(i) We tune J2/J1 to reproduce the neutron-scattering 
result Qb = 0.173 x 27r/5 [S14, S15, SIS, S19]. 

(ii) As shown for the present case of = 0.173 x 27r/6 
in Fig. S3, there exists a single curve in the space 
of {D/Jj^,e,D'JJj_) that reproduces {ipi,ip2)/{2'K) = 
(—0.08,-0.02), which is determined so as to fit the ^Li 
NMR spectra [S15, S20, S21] (see the next section). 
Thus, D/J± and 9 are fixed to satisfy this condition. 

(iii) The rest of the parameters, i.e., Ji, Ja, J±, 
and D'^, are determined to explain the experimental re- 
sults on the anisotropy in the antiferromagnetic reso- 
nance spectra [S16], Im Sx{^) — Im [x'^'^i'^) — (x"°(^) + 
x''''(aj))/2]. The energy levels of the modes observed with 
the THz spectroscopy [S16] are given by Ea — 0.83, 
Eb - 1.04, {Eci,Ec2) = (1.43,1.51), Ed = 2.39, and 
{Eei,Ee2) = (2.98,3.06) in meV. We tune Ji, Ja, J±, 
and Z?" to reproduce the energy levels and the intensities 
of three main peaks observed for Im i5x(a;) at a; = Ea_, 
Eci, and Eei [S16] [Fig. 3(b) of the main body]. 



■20 



DJ.J 
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D/J 
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FIG. S3: The solutions of (D/ J±, D" / J±, 9) for which the 
classical ground state satisfies the condition (tfii, ip2)/ {2n) — 
(-0.08, -0.02) required for fitting the NMR spectra [Fig. S4]. 
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'LI NMR SPECTRA FOR LICU2O2 



[S15, S20, S21], as follows. 



Two parameters (pi {— ipi — —ipm) and (p2 {— 
if 1 1 — —ipiv) of the magnetic ordering pattern are 
fixed to reproduce the ^Li NMR spectra for LiCu202 



In a magnetically long-range ordered phase, the ^Li 
NMR spectra could be dominated by the dipole hyperfine 
field due to the ordered spin moments; 



hdip{ru) 



E 



r-Li) 



(S25) 



with the g factor of electrons, the Bohr magneton /i^, the 
electronic spin Sm,i,j.i at the Cu^"*" position Vm^ij/, and 
the position ru of a Li nuclei. The NMR spectrum at 
the static applied magnetic field H is calculated through 



liH, v) = Ave 6 

''Li 



27r 



-7rLi|/f + /idip(rLi)| , (S26) 



with the gyromagnetic ratio for the ''Li nuclear magnetic 
moment, 77Li. 

Fittings to the "^Li NMR spectra at the static in-plane 
(a5-plane) and out-of-plane magnetic fields are performed 
by changing (^i/(27r) and (/52/(27r) by 1%. Then, sev- 



eral different sets of the rotation angles, ((^1, (/32)/(27r) ^ 
(-0.03,0.04), (-0.08,0.10), and (-0.08,-0.02), reason- 
ably reproduce the experimental results by Gippius et 
al. [S20], as shown for the case of ((^1, (p2)/(27r) = 
(-0.08, -0.02) in Fig. S4. We have not introduced an ad- 
ditional adjustable parameter for the transferred hyper- 
fine coupling, as in the previous analyses [S15, S20]. The 
magnetic patterns of the first and the second sets with 
Lpi ss —Lp2 are not stabilized in our case of 13" 7^ 0. Thus, 
we take the last set (</?!, ^2)/(27r) = (-0.08,-0.02). 



LINEAR SPIN- WAVE THEORY FOR LICU2O2 



0: = -().(t2 




FIG. S4: The ''Li NMR spectra for the afe-twinned LiCu202 
as functions of the dc magnetic field H at the ac mag- 
netic field with the frequency v — 33.4 MHz for (lys 1,(^92) = 
(-0.08,-0.02). (a) The spectrum averaged over the H \\ a 
and H \\ b configurations, (b) The spectrum for H \\ c. 



Following the conventional linear spin-wave approxi- 
mation for the helical magnetism [S22], we introduce 
the Holstein-Primakoff bosons for the spin-S* operators 
= {T'^,Tv,T^) at the site r as 



1 - ^,bX 



= S-blbr. 



br, (S27a) 

_ 1 

' , (S27b) 
(S27c) 



We rotate Tj. so that its z component represents the pre- 
dicted magnetic structure and its x component is perpen- 
dicular to the spiral plane. Thus, we perform the unitary 
transformation Sr — UrTr to obtain 

= COS ipm - 7"™,, siu (p„, sin 6l„i,,; 

+ , i J- sin COS , i J- , ( S28a) 

(S28b) 

S^m,u,t = - j,f sin - , cos Lp,n siu 6'„,, J- ^ 



, j COS ipm cos Om,i.j,e.- 



(S28c) 



We substitute Eqs. (S27) into Eqs. (S28) and then 
Eqs. (S28) to the Hamiltonian Eq. (SI2), and take only 
the bilinear terms in the bosons. 
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Mechanism 


Magnon 


Intra-chain 
(a«;i) 


Inter-chain 


a-axis (aK,ms;o) 


II-III, IV-I (Q„;x) 


I-II, III-IV (a'^) 


Chirality 


fc = 


NA 


NA 


-E" II a, b, c 


E" II c 


k = ±Q 


1! a,c 


NA 


-E" II o, b, c 


E" II b, c 


Magnetostriction 

(CKms;^ ) 


k = 


NA 


NA 


NA 


NA 


k = ±Q 


NA 


NA 


NA 


E'^ II a, c 



TABLE SI: The selection rules for various contributions from bonds to the electric-dipole {E'^) excitation through the 
chirality and magnetostriction mechanisms. The rows for "Intra-chain", "a-axis", "II-III, IV-I", and "I-II, III-IV" represent the 
contributions from the bonds corresponding to Ji (also J2), Ja, J±, and J' (also D^), respectively. They can have independent 
magnetoelectric couplings are,ms;M- "NA" means no activity, i.e., aK,ms;jj = 0. 



In the Fourier-space representation of the spin-wave 
Hamiltonian, there exist bilinear mixing terms of bosons 
bm,q and bl^^ with the wavevectors q and q±Q. They 
generate higher harmonics in the magnetic structure 
and activate the otherwise silent modes for the optical 
electric-dipole transitions. However, the mixing terms 
arc proportional to inter-layer couplings or DM cou- 
plings, and hence are small compared with the intra- 
layer terms. Therefore, the modification of the spin-wave 
spectrum becomes important only when the energy level 
of the 2D spin-wave dispersion almost vanishes, namely, 
around q = and ±Q [circles in Fig. 3(a) of the main 
body]. This allows us to treat these mixing terms along 
the concept of the degenerate perturbation theory of q 
and q± Q modes for q w 0. Namely, we take the basis 
of 



tostriction to the electric dipole moment [S23, S24] , 



bm,q+Q, _q_Q}m=/,- ,iv), 



(S29) 



and express the linearized Hamiltonian as 'i'^i/g^^^'g 
with q « 0. Then, we perform the Bogoliubov transfor- 
mation by solving the generalized eigenvalue problem, 



rrLSW r> D I? 

Lriiq iSq — rSqUi, 

BlGBq = G, 



(S30) 
(S31) 



where G and E are the 24 x 24 diagonal matrices whose 

elements are given by (1,— 1,1,— I,---) and pairs of lev- 
els, (wi, — wi, a;2 5 — W2, • • • )) respectively. Then, we can 
describe the spin-wave eigenmodes in terms of the origi- 
nal spins by using the inverse Bogoliubov transformation 
and the Fourier transform of the inverse of the uni- 
tary transformation Ur- 



MAGNETOELECTRIC COUPLINGS AND THE 
OPTICAL SELECTION RULES FOR LICU2O2 

Here, we summarize the selection rules for local mag- 
netoelectric couplings in various bonds through both con- 
tributions from the vector spin chirality and the magne- 



r — r 



r — r 



- X {SrX Sr') 



PL = J2''L{ry){Sr-Sr.'), 



,(S32) 
(S33) 



respectively, with i being the index for the spatial di- 
rection. The coefficients ajj(r,r') and a^s(r,r') depend 
on the details of both electronic and ionic polarizability, 
and quantitative microscopic calculations of these coeffi- 
cients are beyond our scope. However, the contributions 
from various bonds can be analyzed from the symmetry 
consideration. The results arc summarized in Table SI. 

As far as the one-magnon contributions are concerned, 
the magnetostriction mechanism is active only for the 
magnons at the q = ±Q but not at the Brillouin zone 
boundary in the case of LiCu202, because the primitive 
unit cell contains only one Cu^"*" ion in both the a and 
b directions. Then, this mechanism works only in the 
bonds connecting the layers I and II and those connecting 
HI and IV, having a negligibly small exchange coupling 
J' between the nearest-neighbor spins. Therefore, the 
magnetoelectric coupling through the exchange-striction 
on these bonds should also be negligibly small. 

In contrast, the chirality contributions arc generally 
effective except the bonds along the a axis. In partic- 
ular, contributions from the inter-layer bonds between 
the layers II and III and between I and IV are geomet- 
rically advantageous from the viewpoint of the quantum 
chemistry. In the first order in the relativistic spin-orbit 
coupling, it can flip the spin at the Cu^+ site in one chain, 
and simultaneously change the orbital shape from dx^_y2 
to dyz or dzx- These t2g orbitals form dp-K bonds with 
the 02p orbital located along the c axis, thus allowing for 
a virtual electron transfer between them. This 02p or- 
bital also has a dpa hybridization with the dx2_y2 orbital 
at the Cu^+ site in the other chain. In this way, the in- 
verse Dzyaloshinskii-Moriya mechanism works efficiently. 
Therefore we took into account only these chirality con- 
tributions and adjusted the coefficients in Eq. (S32) in 
order to reasonably reproduce the experimental data on 
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the anisotropy in the dielectric function obtained from 
the THz spectroscopy [S16]. 
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